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by 
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ABSTRACT 

V 

In  this  report,  an  effective  procedure  for  effecting  on  ARMA 
spectral  model  of  a  time  series  is  described.  This  procedure  is 
predicated  on  minimizing  as  set  of  "basic  error"  terms  as  generated 
from  an  ARMA  model  that  is  hypothesized  as  characterizing  the  time 
series  under  analysis.  The  spectral  estimation  performance  achieved 
in  using  this  approach  has  been  empirically  found  to  be  generally 


superior  to  that  obtained  using  such  contemporary  methods  as  maximum 
entropy  and  the  periodogram. 
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I.  INTRODUCTION 


In  this  report,  a  description  of  the  high  performance  ARMA  spectral 

estimator,  as  in-part  developed  under  the  AFOSR  contract  80-0069,  shall 
be  given.  To  begin,  the  main  objective  of  the  project  was  that  of  deve¬ 
loping  an  effective  method  for  estimating  the  ak  and  bfc  coefficients 
governing  the  ARMA  model 

P  Q 

x(n)  +  l  akx(n  -  k)  =  l  bke(n  -  k)  (1) 

k=l  k-0 

in  which  the  excitation  time  series  (e(n)}  is  taken  to  be  white.  These 
coefficients  are  to  be  selected  so  that  this  ARMA  model  is  "most  consis¬ 
tent"  with  a  set  of  time  series  observations. 

x(l ) ,x(2) ,* ’ * ,x(N) .  (2) 

which  are  available.  The  term  "most  consistent"  is  here  being  used  in  the 
sense  that  the  hypothesized  ARMA  model  is  most  compatable  with  the  task  of 
minimizing  a  set  of  "basic  error  terms".  A  brief  description  of  this  pro¬ 
cedure  shall  be  now  given  while  a  more  detailed  description  is  to  be  found 
in  the  appendix. 


I I .  DESCRIPTION  OF  ARMA  MODEL I NG  METHOD 


The  basis  for  the  ARMA  method  Is  dependent  on  the  so-called  basic 
error  terms.  These  terms  are  generated  by  first  multi plyinq  each  side 
of  relationship  (1)  by  the  delayed  entity  x*(n  -  M)  to  obtain 

P 

e(m,n)  =  x(n)x*(n  -  m)  +  l  akx(n  -  k)x*(n  -  m)  (3) 

k=l 

q  +  1  £  m  <  N 
max  (m,p)  £  n  £  N 

It  can  be  shown  that  these  basic  error  terms  e(m,n)  are  each  zero  mean 
when  the  underlying  time  series  in  an  ARMA  process  of  order  less  than  or 
equal  to  (q,p)  and  the  ak  coefficients  in  expression  (3)  correspond  to 
the  exact  process',  autoregressive  coefficients. 

With  this  in  mind,  a  method  for  selecting  the  ARMA  models  ak 
coefficients  is  suggested.  Namely,  they  are  chosen  so  as  to  cause  the 
basic  error  terms  (3)  to  be  as  close  to  their  mean  value  of  zero  as 
possible.  This  objective  may  be  realized  by  introducing  the  following 
quadratic  functional 

f(a)  =  e-tWe  (4) 

in  which  £  is  a  column  vector  whose  elements  are  appropriate  arrangements 
of  the  ensemble  of  basic  error  terms  (3),  W  is  a  nonegative  matrix,  a  is 
the  autoregressive  coefficient  vector  with  elements  ak,  and,  +  denotes  the 
operation  of  complex  conjugate  transposition.  It  is  readily  shown  that 
the  minimization  of  functional  (4)  will  result  in  a  set  of  consistent 


linear  equations  for  the  optimum  ARMA  model  autoregressive  coefficients. 

Once  the  afc  coefficients  have  been  obtained,  the  ARMA  model's 
moving  average  coefficients  are  estimated  by  first  generating  the  so- 
called  residual  sequence  (e(n)}  according  to 

P 

e(n)  =  x(n)  +  £  akx(n  "  p  +  1  <  n  <  N  (5) 

k=l 

The  spectrum  of  this  residual  time  series  is  theoretically  given  by 

q 

|  l  bke"'I*cw  |^.  A  technique  for  obtaining  a  moving  average  estimate 
k=0 

of  this  residual  time  series  is  fully  described  in  the  appendix  and 
entails  the  utilization  of  the  smoothed  periodogram.  Once  this  smoothed 
periodogram  has  been  generated,  the  required  ARMA  spectral  estimate  is 
achieved.  Empirically  derived  results  indicate  that  this  report's  pro- 
cedure  provides  a  superior  spectral  estimation  performance  than  that 
achieved  by  such  contemporary  alternatives  as  the  maximum  entropy  method, 
and  the  periodogram. 


III.  CONTRACT  PUBLICATIONS 

The  following  two  refereed  publications  resulted  from  the  sponsored 
AFOSR  contract. 

[1]  "  Autoregressive  Moving  Average  Spectral  Estimation:  A  Model 
Equation  Error  Procedure",  IEEE  Trans,  on  Geoscience  and 
Remote  Sensing,  vol.  GE-19,  No.  1,  Jan.  1981,  pp  24  -  28. 

[2]  "Two  Dimensional  Spectral  Estimation",  IEEE  Trans,  on  Acoustics, 
Speech  and  Signal  Processing,  vol.  ASSP-29,  No.  3,  June  1981, 

pp  396  -  401 . 
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Autoregressive  Moving  Average  Spectral  Estimation: 
A  Model  Equation  Error  Procedure 

JAMES  A.  CADZOW,  senior  member,  ieee 


Abtawct-k  procedure  it  pramnwd  for  growl  ling  an  lutongrwive 
movag  magi  (ARMA)  ipactral  modal  of  a  stationery  lima  saties 
haaad  apoa  a  Borta  sat  of  tfana  aariaa*  obaarratioar,  Tht  ARMA  modal's 
aatoiapamhn  coafflciaats  are  aaUmatad  by  minim  Bing  a  quadratic 
ftmctfcm  of  a  aat  of  baac  error  terms.  In  exam  pies  treated  to  data,  this 
method  haa  demonstrated  an  exceptional  ability  in  re  solving  closely 
spaced  narrow  band  atomic  in  a  low  dpaai-to  nows  environment  where 
other  procadntaa  each  aa  the  maximum  entropy  method  often  fail.  Its 
effectiveness  on  othar  dames  of  time  series  also  dtows  promim  and  a 
more  pnand  evdaetkin  to  presently  being  conducted.  With  this  in 
mind,  the  now  ARMA  procedure  prom  tom  to  be  an  important  ro*ctral 
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I.  Introduction 

N  this  paper,  we  shall  be  concerned  with  the  task  of  esti¬ 
mating  the  statistical  characteristics  of  a  stationary  random 
time  series  {x(n)}  from  a  finite  set  of  observations.  For  many 
applications,  knowledge  of  the  time  series’  underlying  autocor¬ 
relation  sequence  as  formally  defined  by 

r,  («)«£{*(;.  ♦*)*•(*)}  (1) 

conveys  all  the  information  required.  In  this  expression,  £  and 
*  denote  the  expected  value  and  complex  conjugation  opera¬ 
tions,  respectively.  It  is  often  advantageous  to  equivalently 
characterize  stationary  time  series  in  the  frequency  domain 
where  their  intrinsic  properties  may  be  more  discernible.  This 
is  particularly  true  for  so-called  narrow-band  processes.  The 
vehicle  for  this  characterization  is  the  associated  power  spec- 
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tral  density  as  given  by 

S,(w)=  t  rx(n)e-^n  (2) 

«■-* 

which  is  recognized  as  being  the  Fourier  transform  of  the 
autocorrelation  sequence. 

Upon  examination  of  expressions  (1)  and  (2),  it  is  apparent 
that  the  required  time  series  characterization  will  generally  en- 
tail  complete  knowledge  of  the  generally  “infinite”  length 
autocorrelation  sequence.  Here,  we  shall  be  concerned  with 
extracting  this  information  from  the  following  “finite”  set  of 
time  series  observations 

x(l),*(2),  •••,*(#).  1  (3) 

Unless  some  constraints  are  imposed  upon  the  time  series' 
basic  nature,  it  is  clear  that  there  exists  a  fundamental  incom- 
patability  in  estimating  the  required  statistical  knowledge  from 
this  finite  set  of  measurements.  In  accordance  with  well  ac¬ 
cepted  practices,  this  dilemma  is  here  resolved  by  postulating  a 
finite  parameter  linear  model  to  represent  the  random  phe¬ 
nomenon.  Specifically,  the  response  of  this  model  to  “ran¬ 
dom”  shock  excitations  will  be  conceived  of  as  generating  the 
time  series  under  study. 

In  terms  of  parameter  parsimony,  the  causal  autoregressive 
moving  average  (ARMA)  model  of  order  (p,q)  as  specified  by 

x(n)  +  £  akx(n  -  k)  =  £  bke(n  -  k)  (4) 

k - 1  k -0 

is  generally  the  most  effective  linear  model.  In  this  representa¬ 
tion,  the  unobserved  excitation  (e(/i)}  is  taken  to  be  a  white 
noise  time  series  with  zero  mean  and  variance  a2.  It  is  impor¬ 
tant  to  appreciate  the  fact  that  the  more  specialized  autore¬ 
gressive  (i.e.,  bk  =  0  for  k  0)  linear  model  will  generally  en¬ 
tail  a  much  higher  order  choice  so  as  to  achieve  a  comparable 
statistical  representation.  Conceptually,  the  more  efficient 
ARMA  model  is  the  logical  model  choice  when  the  exact  na¬ 
ture  of  the  time  series  is  unknown. 

Due  to  the  relatively  difficult  task  of  estimating  the  ARMA 
model’s  ak  and  bk  parameters  from  the  given  observations  (3), 
however,  the  preponderance  of  activity  has  been  directed  to¬ 
wards  the  specialized  autoregressive  (AR)  model.  This  is  a  di¬ 
rect  consequence  of  the  simpler  AR  parameter  estimation 
problem  which  results.  In  particular,  the  basically  equivalent 
maximum  entropy,  linear  predictive  coding,  and  AR  methods 
have  been  developed  for  efficiently  obtaining  the  AR  model 
parameters  [  1  ] .  Nonetheless,  the  inherent  superiority  of  an 
ARMA  model  representation  is  widely  recognized  and  a  num¬ 
ber  of  procedures  for  estimating  the  ARMA  model  parameters 
have  been  advanced.  These  include  the  whitening  filter  ap¬ 
proach  which  is  iterative  in  nature,  generally  slow  in  conver¬ 
gence,  and  typically  requires  a  relatively  long  data  length  N  to 
be  effective  (2)  and  [3] .  A  more  desirable  closed-form  ap¬ 
proach  which  is  based  on  the  autocorrelation  relationship 
governing  ARMA  processes  has  also  been  developed  [4]  -[6] . 
This  approach  often  provides  good  estimates  and  does  not  re¬ 
quire  jV  to  be  excessively  large.  Unfortunately,  its  perfor¬ 
mance  is  not  always  as  good  as  that  provided  by  AR  methods. 


In  this  paper,  a  generalization  of  a  recently  developed  closed- 
form  ARMA  method  shall  be  presented  [7]  and  [8].  This 
new  procedure  has  been  empirically  found  to  provide  signifi¬ 
cantly  better  spectral  estimation  performance  than  the  above 
two  ARMA  approaches  as  well  as  the  maximum  entropy 
method.  In  what  is  to  follow,  a  time  domain. approach  for  de¬ 
termining  the  ARMA  model's  ARa*  coefficients  will  be  given. 
This  in  turn  will  be  followed  by  a  frequency  domain  proce¬ 
dure  for  estimating  the  effects  of  the  moving  average  bk  coef¬ 
ficients  on  the  overall  spectral  estimate.  Use  will  be  made  of 
the  well-known  fact  that  the  power  spectral  density  corre¬ 
sponding  to  ARMA  model  (4)  is  specified  by 


|1  +ate-‘"  +  --  +  ape-*'"\i 


O' 


(5) 


II.  Autoregressive  Coefficient  Estimates 
An  effective  method  for  estimating  the  ARMA  model's  AR 
coefficients  from  the  set  of  given  observations  will  be  now  pre¬ 
sented.  This  first  entails  multiplying  each  side  of  expression 
(4)  by  the  entity  x*(n  -  m )  to  yield  the  “basic  error  terms” 


p  : 

e(m,  n)  *  x(n)  x*(n  -  m)  +  £  <***(«  '  k)  x*(n  -  m)  (6a) 

k  •  1 
<7 

=  £  bke(n  -  k)x*(n  -  m),  .  for  q+l<m<;V 
*«o 

max  (i m,p)<n<N 
(6b) 

where  the  range  on  the  m  and  n  variables  is  dictated  by  the 
given  time  series  range  of  1  <n  <N  and  the  desire  to  have  the 
random  variables  e(n  -  k)  and  x*(n  -  m)  be  uncorrelated.  If 
the  time  series  is  in  fact  an  ARMA  process  of  order  less  than  or 
equal  to  (p,q)  it  follows  that  the  basic  error  terms  e(m,  n)  are 
each  zero  mean  random  variables.1  This  is  a  consequence  of 
the  causality  of  model  (4)  and  the  whiteness  of  the  excitation 
which  results  in  the  random  variables  e(n  -  k)  and  x*(n  -  m) 
being  uncorrelated. 

With  these  thoughts  in  mind,  relationship  (6)  provides  an 
ideal  vehicle  for  determining  a  set  of  AR  coefficients  which  are 
consistent  with  the  given  time  series  observations.  Specifi¬ 
cally,  the  ak  coefficients  will  be  selected  so  as  to  cause  each  of 
the  basic  error  terms  (6a)  to  be  as  close  to  their  mean  value  of 
zero  as  possible.  This  goal  is  achieved  by  minimizing  a  squared 
error  criterion  of  the  form 

/(«)  =  e+We  (7) 

in  which  e  is  a  column  vector  whose  elements  are  appropriate 
arrangements  of  the  ensemble  of  basic  error  terms  f6a).  tV  is  a 


1  As  a  side  note,  it  may  be  shown  that  the  basic  error  terms  have  iden¬ 
tical  variances  given  by 

1 

fjc(O)  a2  £  |h*lJ. 

Ac*  0 
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nonnegative  definite  square  matrix,  and,  t  denotes  the  opera¬ 
tion  of  complex  conjugate  transposition.  This  criterion  is 
readily  shown  to  be  a  function  of  the  AR  coefficients  upon 
substitution  of  expression  (6a)  into  criterion  (7). 

The  weighting  matrix  W  is  to  be  selected  based  on  statistical 
reasoning  which  must  be  tempered  by  an  appreciation  that  lit¬ 
tle  is  generally  known  about  the  time  series’  statistics.  To  illus¬ 
trate  this  point,  let  us  consider  the  usually  hypothesized  condi¬ 
tion  that  the  excitation  {e(n)}  is  a  Gaussian  white  process.  In 
this  case,  if  the  matrix  W  is  set  equal  to  E  {eet}‘1 ,  the  minimi¬ 
zation  of  criterion  (7)  will  result  in  a  maximum  likelihood  esti¬ 
mate  of  the  AR  coefficients.  Unfortunately,  the  generation  of 
this  particular  matrix  requires  knowledge  of  the  time  series’ 
{x(n)}  second-order  statistics.  This  statistical  information, 
however,  is  precisely  what  we  are  seeking  to  estimate  via  the 
ARMA  model  (4). 2  This  dilemma  has  arisen  due  to  the  unreal¬ 
istic  desire  to  seek  a  maximum  likelihood  ARMA  model. 
More  modest  objectives  must  be  sought  if  a  tractable  proce¬ 
dure  is  to  be  evolved. 

A  reasonable  implementation  of  the  above  autoregressive 
coefficient  selection  philosophy  then  depends  on  a  prudent 
choice  for  the  weighting  matrix.  In  this  paper,  we  will  be  con¬ 
cerned  exclusively  with  a  choice  that  results  in  the  following 
mean-squared  error  criterion: 


f(a)  =  £  w(m) 


m  »q  + 1 


N  2 

Z  e(m,  n) 

n  -  l  +m» x(m,p) 


(8) 


although  other  choices  are  suggested  in  Section  V.  In  this  ex¬ 
pression,  the  w(m)  are  nonnegative  parameters  which  are  logi¬ 
cally  selected  to  be  inversely  proportional  to  the  variance  of 
the  term  which  they  multiply.  A  general  expression  for  these 
variance  entities,  unfortunately,  also  reveals  a  dependency  on 
the  time  series’  unknown  second-order  statistics  and  the 
ARMA  model  coefficients.  Empirical  evidence  suggests,  how¬ 
ever,  that  a  choice  of  the  weighting  parameters  as  given  by 

w(m)  =  [N  -  m]*  q  +  l<m<N-l  (9) 


results  in  satisfactory  performance.  An  examination  of  other 
weighting  parameter  selections  is  currently  being  conducted. 

Using  standard  calculus,  it  is  readily  shown  that  the  set  of 
AR  coefficients  which  render  criterion  (8)  a  minimum  must 
satisfy  the  following  linear  system  of  equations: 


Co  *c  (10) 

where  a  is  the  p  X  1  AR  coefficient  vector  estimate  with  ele¬ 
ments  a*,  C  is  a  p  X  p  nonnegative  definite  Hermitian  matrix 
with  elements 


N- 1  N  N 

c«  *  Z  Z  Z  **'(«)  x(»  -  *)*(*  - m) 

m  •q  *  i  /i«i  )■( 

•  x*(n  -  m)x*(s  -  /'),  for  i,k  *  1, 2,  •  •  •  ,p  (11a) 


:!n  the  special  case  of  an  AR  model,  however,  it  is  possible  to  obtain 
the  maximum  likelihood  estimate. 


where  t  *  1  +  max  (m,p)  and  c  is  a  p  X  1  vector  with  elements 
c,*-cl0.  for  l<i<p.  (lib) 

Thus  the  optimum  set  of  ARMA  autoregressive  coefficient 
estimates  are  obtained  by  solving  the  linear  system  of  equa¬ 
tion  (10).  It  is  possible  to  utilize  the  Hermitian  structure  of 
matrix  C  to  yield  an  efficient  solution  procedure. 

A.  Forward  and  Backward  Data  Approach 
One  can  more  effectively  use  the  observed  data  (3)  to  achieve 
an  ARMA  spectral  estimate  by  using  that  data’s  backward  ver¬ 
sion  as  specified  by 

x(1).*(2),  -.*(A0  (12) 

in  which  x(n)sx*(N  +  1  -  n).  In  particular,  one  formulates  a 
mean-squared  error  criterion  as  given  by 

N-l  f|  N  2 

/(*)=  Z  Z  e(m,n) 

m  ~q  *  i  t|n  ■  l  ♦  max  (m.p) 


N  2 1 

Z  ?(m* ")  r 

•  l  ♦  max  (m.p)  J 


(13) 


where  e(m,  n )  denotes  the  basic  error  term  associated  with  the 
backward  sequence  (12).  It  is  a  simple  matter  to  show  that 
the  AR  coefficient  vector  which  minimizes  this  forward- 
backward  criterion  must  satisfy  the  linear  system  of  equations 

[C  +  C]  a  -  c+c  (14) 

where  C  and  C  are  each  pXp  nonnegative  definite  Hermitian 
matrices  given  by  relationship  (11a)  with  the  forward  (3)  and 
backward  (12)  data  used,  respectively,  and,  the  p  X  1  vectors 
c  and  c  are  given  similar  interpretations. 

It  is  important  to  note  that  the  rank  of  matrix  [C  +  C]  at 
least  equals  the  rank  of  either  of  its  constituent  nonnegative 
definite  matrices  C  or  C.  Thus,  in  using  this  forward  and  back¬ 
ward  approach  to  estimate  the  AR  coefficients,  one  is  typi¬ 
cally  able  to  use  a  higher  order  model  (i.e.,  larger  value  of  p) 
than  would  be  the  case  for  the  forward  or  backward  models 
only.  Moreover,  it  has  been  empirically  found  that  the  for¬ 
ward  and  backward  approach  usually  results  in  better  spectral 
estimates. 


III.  Numerator  Dynamics 

A  variety  of  procedures  exist  for  determining  the  numerator 
dynamics  of  an  ARMA  time  series  once  the  AR  coefficients 
have  been  estimated  {4 J  -[9] .  In  this  section,  a  procedure 
which  has  been  found  to  produce  improved  results  will  be  pre¬ 
sented.  It  makes  use  of  the  characteristic  equation  (4)  which 
can  be  interpreted  as  generating  the  auxiliary  “residual’’  se¬ 
quence  (e(n)}  according  to 

p  c 

e(n)sx(n)  +  £  akx(n-k)  '  (15a) 

*■1 

=  Z  bke(n-k).  (15b) 

*-o 

One  may  straightforwardly  show  that  the  spectrum  of  the  mov¬ 
ing  average  residual  time  series  Se(u>)  is  given  by  |R(w)|2o2. 
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This  observation  in  conjunction  with  relationship  (5)  provides 
the  vehicle  for  obtaining  the  underlying  time  series  spectral  es¬ 
timate,  that  is 

S,(w)  =  Se(w)/|^M|l.  (16) 

With  this  in  mind,  a  method  will  be  now  given  for  obtaining  a 
<7th-order  moving  average  spectral  estimate  of  the  residual  se¬ 
quence.  This  estimate  will  be  based  on  the  specific  residual 
elements 

e(p+l),e(p  +  2),--,e(A0  (17) 

which  are  calculated  using  relationship  (15a)  and  the  given 
time  series  observations  (3).  The  AR  coefficient  used  in  (ISa) 
will  correspond  to  the  ak  estimate  obtained  upon  solving  ex¬ 
pression  (14). 

The  approach  to  be  now  presented  in  an  adaption  of  the 
well-known  method  of  Welch  for  obtaining  smoothed  peri- 
odogram  estimates  [10].  In  essence,  one  first  segments  the 
calculated  residual  elements  (17)  into  L  segments  each  of 
length  q  +  1  as  specified  by 

ek(n)  -  a(n)  e(n  +  1  +p  +  fcf),  0  <n<q 

0<k<L  -  1  (18) 

where  a(n)  is  a  data  window  and  “d"  a  positive  integer  which 
specifies  the  time  shift  between  adjacent  segments.  These  in¬ 
dividual  segments  will  overlap  for  a  shift  selection  of  d<q. 
Furthermore,  to  include  only  observed  data,  the  relevant  pa¬ 
rameters  must  be  selected  so  that  p+q  +  l+(L-  1 )  d  <N. 
Finally,  the  periodogram  of  each  of  the  L  segments  is  taken 
and  these  are  averaged  to  obtain  the  desired  smoothed  qih- 
order  periodogram,  that  is 

~  i  i-i  f  i  4  ,1*1 

7  £  £  *(*)*(* +p  +  *tf)-'H 

L  km  ft  I<7+1  nTft 


The  data  window  is  to  be  normalized  according  to  £  ar*(n)  =  1 . 
In  using  this  segmentation,  the  variance  of  the  estimate  5«(<o) 
is  reduced.  The  price  paid  for  this  reduction,  however,  is  a  loss 
in  frequency  resolution,  an  increased  bias  of  the  estimate,  and, 
a  possible  deterioration  in  power  level  estimates.  It  must  be 
noted,  however,  that  the  overall  resolution  capability  of  this 
ARMA  procedure  is  predominately  influenced  by  the  AR  co¬ 
efficient  selection.  If  one  is  basically  interested  in  resolution 
performance,  an  examination  of  the  ARMA  model’s  pole  loca¬ 
tions  then  need  only  be  considered. 

IV.  Numerical  Example 

In  order  to  compare  the  effectiveness  of  the  new  ARMA 
spectral  estimator  with  the  maximum  entropy  method,  the" 
classical  problem  of  resolving  two  closely  spaced  (in  fre¬ 
quency)  sinusoids  in  white  noise  will  be  considered.  Specifi¬ 
cally,  the  time  series  under  study  is  specified  by 

jc(n)  *  y/20  cos  (0.4rrn)  +  s/l  cos  (0.426rrn)  +  w(n)  (20) 

where  {w(n)}  is  a  white  Gaussian  noise  process  of  zero  mean 
and  unity  variance.  The  sinusoids  of  normalized  frequencies 
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Fig.  1.  Fifteenth-order  ARMA  spectral  estimate  of  two  sinusoids  at 
normalized  frequencies  f,  *  0.4  (10  dB)  and  f}  »  0.426  (0  dB)  baaed 
upon  individual  sequences  of  length  64. 
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Fig.  2.  Thirtieth-order  maximum  entropy  AR  spectral  estimate  of  two 
sinusoids  at  normalized  frequencies/]  »  0.4  (10  dB)  and  /j  ■  0.426 
(0  dB)  based  upon  individual  sequences  of  length  64. 

0.4  and  0.426  are  readily  calculated  to  have  a  signal-to-noise 
ratio  (SNR)  of  10  dB  and  0  dB,  respectively.  A  sequence  of 
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Fig.  3.  ARMA  spectral  estimates  of  two  sinusoids  at  normalized  fre¬ 
quencies  l\  ■  0.4  (10  dB)  and  *  0.426  (0  dB)  based  upon  all  640 
observations. 
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Fig.  4.  Thirtieth-order  maximum  entropy  AR  spectral  estimate  of  two 
sinusoids  at  normalized  frequencies  f\  ■  0.4  (10  dB)  and  /j  *  0.426 
(0  dB)  based  upon  all  640  observations. 

length  640  defined  over  0  <n<  639  was  next  generated  using 
this  relationship.  Furthermore,  in  order  to  provide  a  statistical 
basis  for  our  comparison,  this  640  length  sequence  was  then 
decomposed  into  ten  disjoint  sequences  each  of  length  64  de¬ 
fined  on  0  </i  <63,  64<n  <  127,  •  •  • ,  576  <n  <639.  An 
ensemble  consisting  of  ten  subsequences  each  of  length  64  has 
thereby  been  generated  with  each  subsequence  having  a  differ¬ 
ent  noise  sample  and  a  different  initial  phase  between  the  two 
sinusoids.  This  latter  condition  is  useful  in  revealing  any  po¬ 
tential  sensitivity  to  initial  phase  that  the  new  ARMA  spectral 
estimation  method  may  possess. 

The  spectral  estimates  which  resulted  when  the  new  ARMA 
spectral  estimator  (forward  data  only)  and  the  maximum  en¬ 
tropy  (covariance)  method  were  applied  to  these  ten  random 
time  series  samples  are  displayed  in  Figs.  1  and  2,  respectively. 
The  ordinates  were  scaled  from  0  to  50  dB  for  each  of  the  ten 
individual  plots  with  each  AR  maximum  entropy  estimate  be¬ 
ing  of  order  thirty  and  each  ARMA  estimate  being  of  order  fif¬ 
teen  (with  d*  16).  The  spectral  estimates  are  plotted  one 
above  the  other  (lower  one  corresponds  to  0  <  n  <  63)  so  as 
to  better  reveal  the  consistency  of  the  estimates  and  to  con¬ 
serve  space.  It  is  clear  that  the  maximum  entropy  method  is 
able  to  resolve  the  sinusoids  in  only  five  of  the  ten  sample 
cases.  On  the  other  hand,  the  new  ARMA  method  successfully 
resolves  the  sinusoids  in  each  case  and  provides  excellent  fre¬ 
quency  estimates  of  fx  »  0.3980  and  /2  *  0.4295  obtained  by 
averaging  over  the  ten  cases  with  associated  sampled  standard 
deviations  of  a/(  =  0.0027  and  af%  *  0.0107. 

To  illustrate  the  new  ARMA  method’s  effectiveness  on  long 
data  length  sequences,  an  eighth-  and  fifteenth-order  ARMA 
model  were  next  obtained  using  the  entire  sequence  of  length 
640  as  given  above.  The  spectral  estimates  which  resulted  are 
displayed  in  Fig.  3  where  a  sharpness  in  frequency  resolution  is 


apparent.  In  addition,  a  thirtieth-order  AR  spectral  estimate 
which  arises  from  this  640-length  sequence  is  shown  in  Fig.  4. 
Gearly,  the  new  ARMA  method  has  outperformed  the  AR 
maximum  entropy  method  in  both  the  short  and  long  data 
length  examples  here  considered. 

V.  Conclusions 

A  generalization  of  a  recently  developed  ARMA  spectral  esti¬ 
mation  method  has  been  presented.  This  has  included  the  in¬ 
troduction  of  an  error  weighting  matrix,  the  concept  of  using 
forward  and  backward  data,  and,  the  utilization  of  Welch’s 
method  for  obtaining  estimates  of  the  spectrum’s  numerator 
dynamics.  Empirical  evidence  suggests  that  this  new  proce¬ 
dure  has  an  improved  spectral  resolution  performance  when 
compared  to  popular  contemporary  methods.  Its  full  poten¬ 
tial  will  be  realized,  however,  only  after  a  number  of  relevant 
issues  are  resolved.  These  include  the  task  of  determining 
good  choices  for  the  weighting  matrix  W,  obtaining  a  data  de¬ 
pendent  procedure  for  estimating  the  ARMA  model  order, 
and,  investigating  other  numerator  dynamics  methods. 

In  addition  to  the  specific  weight  matrix  selection  consid¬ 
ered  in  this  paper,  the  following  two  choices  of  criterion 

N-l  N 

£  £  \e(m,n)\2 

m  *q  + 1  m  *  m  ♦  1 

N-l  N-l  2 

and  £  w(n)  ]T  e(m,n) 
nmq* 2  m*q* l 

have  given  preliminary  evidence  of  providing  satisfactory  spec¬ 
tral  estimation  performance.  Further  investigation  of  this 
most  critical  weighting  matrix  choice  is  now  being  conducted 
and  will  be  subsequently  reported  upon. 
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